 function [GG,C,RR,eu,SDX,ZZ,NY,NX,addout]=modchimodDtFacpInfExpSplit_obs(parest,options,op,calval,posest,poscal)
% =========================================================================
%
% function [GG,C,impact,eu,SDX,ZZ,NY,NX,NETA,gev,ncof]=...
% modchimodDtFacpInfExp_obs(parest,calval,posest,poscal)
%
%% *I. Description for CHIMOD DT FACP INFEXP _ OBS * 
% 
% This is the OBS version of the code used for estimation 
% Description for CHIMOD DT FACP INFEXP
% 
% Model code for Chicago Fed model (CHIMOD) with 
% Double Trend (DT) Factor structure for prices (FACP) 
% Separates the ARMA(1,1) wage markup/labor disutility shock into a 
% white noise wage markup and an AR(1) labor disutility shock 
% with a spread observable and allowing for policy news shocks
% upto 4 periods ahead, Incorporates CPI inflation and SPF inflation 
% expectations (INFEXP) observables with inflation target in policy rule . 
%
% Refer to that code for all pertinent documentation. Here will only focus
% on inputs and outputs that are different. 
% 
%% *II. Input* 
%
% There are only 2 inputs 
% *PARAM:*  Complete vector of parameters 
% *ADDIN:*  Structure with additional inputs 
% 
% .STDSIGNALS  [12 1] vector of STD for the MP signals 
%              Any > 0.01 will be treated as active. 
%              If not provided, STD will be set to 1e-5 
% 
%% *III. Output*  
%
% The model to be solved is written, inside this code, as $$ A_{0}Y_{t+1} +
% A_{1}Y_{t}=A_{2}Y_{t-1} + PSI \eta_{t} $$  where the A matrices are
% filled with the linearized equation and entries depend on the parameter
% vector PARAM. 
%
% The solution of the model is a state space system characterized by 
%
% *Transition Equation* $$ S_{t} = GG S_{t-1} + RR \eta_{t} $$ and variance covariance matrix, 
% where $$ S_{t} $$ is the vector of all
% variables (endogenous and exogenous) of the model of dimension [NY 1],
% and $$ \eta $$  is the vector of exogenous innovations of dimension [NX 1] with 
% variance covariance matrix $$ V(\eta) = SDX'SDX $$. 
% Matrix dimensions: 
% GG [NY NY], IMPACT [NY NX], SDX [NX NX].  
%
% *Measurement Equation* $$ Y_{t} = ZZ S_{t} + ZZ C  $$ , 
% where $$ Y_{t} $$ is the data vector of dimension [NZ 1], $$ ZZ $$ is a selection matrix indicating 
% which states are observed, and $$ CC $$ is a vector of constants of
% dimension [NY 1]. 
% 
% EU is a [2 1] vector equal to [1;1] if there is a unique stable RE
% solution 
% 
%% III.a STRUCTURE ADDOUT 
%
% Relative to estimation code there are numerous new outputs, in the 
% structure ADDOUT. 
%
% III.a Additional Matrices 
% 
% *.Z2:*            The model observation equation will be enhanced by the matrix $$ Z_{2} $$
%                   which should match dimension of additional data not
%                   used for estimation but used for filtering 
% *.C2:*            [NY 1] Vector of constants. Includes all constants from
%                   the observation equation, plus additional constants for all states with 
%
% *.STANAMES:*      [NY 1] Cell with the names of all States with a SS
% mean. *Note:* to demean the additional observables can then use $$ Z2*CC2
% $$. 

% *.STSHOCKS:*      [NX 1] Cell with the names in STANAMES that correspond to
%                      shocks 
% *.SHONAMES:*      [NX 1] Cell with the names of the innovations (will be used
%                      in graphs)
% 
%% *Version*
% Created by Scott Brave 07/13/2011  
% 
% =========================================================================

if nargin~=1 && nargin~=2 
    error('Need either 1 or 2 inputs') 
end 

% =========================================================================
%% *PART 1. Declare All Variables in the Model Solution $$ S_{t} $$* 

%__________________________________________________________________________
% 1.a Endogenous Variables 
y=1;            % output
k=2;            % capital
L=3;            % hours
Rk=4;           % return on capital
w=5;            % real wages
x=6;            % marginal utility of labor
p=7;            % inflation
s=8;            % marginal cost
lambda=9;       % multiplier
c=10;           % consumption
R=11;           % interst rate
u=12;           % capital utilization
phi=13;         % multiplier
i=14;           % investment
kbar=15;        % "gross" capital
%__________________________________________________________________________
% 1.b Exogenous driving processes
mp=16;
z =17;          % Productivity shock
g =18;          % Public spending shock
miu=19;         % Relative price (non stationary), this is Growth Rate
lambdap=20;     % Markup shock AR(1) 
psi=21;         % Preference shock AR(1) 
b=22;           % Discount factor shock 
upsilon=23;     % MEI, stationary, investment shock
pitarg=24;      % Inflation Drift
wmarkwn=25;     % White Noise Wage Markup 
%__________________________________________________________________________
% 1.c Observable growth rates (nominal) 
ynomobs=26; 
cnomobs=27;
inomobs=28; 
wnomobs=29; 
%__________________________________________________________________________
% 1.c.a Auxiliary Variables or Observable Prices
yrobs   =30; 
gdp     =31; 
dpinv   =32; 
gdpdef  =33; 
crobs   =34;   % delta(Nominal C - Pc)
irobs   =35;   % delta(Nominal I - Pc)
wrobs   =36;   % delta(Nominal W - Pc)
ygap    =37;   % Output Gap 
exreal1 =38;   % Expected Real Rate next period 
%__________________________________________________________________________
% 1.d Lags
yrobs_1=39;   % Real Observable Output (adjusted for utilization) lag1      
yrobs_2=40;   % Real Observable Output (adjusted for utilization) lag2
yobsan =41;   % Annualize Observable Output 
p_1    =42;   % Lagged C deflator 
p_2    =43;   % Lagged C deflator 
dpma   =44;   % Four quarter inflation  
dagtrend=45;  % First difference of the aggregate trend 
% Make DAGTREND Last Variable here for automatic updating of the next
% position 
%__________________________________________________________________________
% 1.e 1 through 12 step ahead expectations of FFR 
startRlead=dagtrend+1;  % Number to start with CURRENTLY 46
R_lead1=startRlead;  
R_lead2=startRlead+1; 
R_lead3=startRlead+2; 
R_lead4=startRlead+3;
R_lead5=startRlead+4;
R_lead6=startRlead+5;
R_lead7=startRlead+6;
R_lead8=startRlead+7;
R_lead9=startRlead+8;
R_lead10=startRlead+9;
R_lead11=startRlead+10;
R_lead12=startRlead+11;
Rlead_vec=(R_lead1:R_lead12); 
%__________________________________________________________________________
% 1.f Policy signals and Inflation Expectations
startSig=startRlead+12; % Number to start with CURRENTLY 58 
sigR_0=startSig; 
sigR_1=startSig+1; 
sigR_2=startSig+2; 
sigR_3=startSig+3; 
sigR_4=startSig+4; 
sigR_5=startSig+5; 
sigR_6=startSig+6; 
sigR_7=startSig+7; 
sigR_8=startSig+8; 
sigR_9=startSig+9; 
sigR_10=startSig+10; 
sigR_11=startSig+11; 
sigR_12=startSig+12; 

posExpInf = startSig+12+1;
pe_vec = (posExpInf):(posExpInf+39);
expectedPi40 = posExpInf+40;
realRate12q  = posExpInf+41; 
%__________________________________________________________________________
% 1.g Flex Price and Wage Economy * 
startstar=posExpInf+41+1;        % Number to start with CURRENTLY 111
ystar=startstar;              % output
kstar=startstar+1;            % capital
Lstar=startstar+2;            % hours
Rkstar=startstar+3;           % return on capital
wstar=startstar+4;            % real wages
xstar=startstar+5;            % marginal utility of labor
Rstar=startstar+6;            % interst rate
sstar=startstar+7;            % marginal cost
lambdastar=startstar+8;       % multiplier
cstar=startstar+9;            % consumption
ustar=startstar+10;           % capital utilization
phistar=startstar+11;         % multiplier
istar=startstar+12;           % investment
kbarstar=startstar+13;        % "gross" capital
gdpstar=startstar+14;         % GDP 
NY=startstar+14;
% =========================================================================

% =========================================================================
%% *PART 2. Declare Innovations to the exogenous variables $$ \eta_{t} $$* 

nsignals = 12;    % Number of policy signals 
NX=10+ nsignals; 
%__________________________________________________________________________
% 2.a Structural Disturbances 
Rs=1;           % MP
zs=2;           % Technology 
gs=3;           % Public spending
mius=4;         % Investment     
lambdaps=5;     % Mark-up
psis=6;         % Leisure preference
bs=7;           % Discount factor
upsilons=8;     % investment shock
shpitarg=9;     % Inflation Drify
wmarkwnsh=10; 
%__________________________________________________________________________
% 2.b Policy Signals 
Rs_lag1=11;    % 1 period ahead
Rs_lag2=12;    % 2 periods ahead
Rs_lag3=13;    % 3 periods ahead
Rs_lag4=14;    % 4 periods ahead
Rs_lag5=15;    % 5 periods ahead
Rs_lag6=16;    % 6 periods ahead
Rs_lag7=17;    % 7 periods ahead
Rs_lag8=18;    % 8 periods ahead
Rs_lag9=19;    % 9 periods ahead
Rs_lag10=20;   % 10 periods ahead
Rs_lag11=21;   % 11 periods ahead
Rs_lag12=22;   % 12 periods ahead
% =========================================================================

% =========================================================================
%% *PART 3. Declare Parameters* 
%__________________________________________________________________________
% 3.1 Parameters for all but multiple prices, GDP and I deflators 
ncof=        2*68;    % Number of coefficients 
numpar=      ncof;  % Equal to Parameters, distinction irrelevant here 
if nargin < 4 || isempty( calval ) 
    param=parest  ; 
else
    param=zeros(numpar,1); 
    param(posest)=parest; 
    param(poscal)=calval; 
end; 

alpha=param(1);         % share of capital in the prod. function
delta=param(2);         % capital depreciation rate
iotap=param(3);         % price indexation (=0 is static indexation, =1 is dynamic)
iotaw=param(4);         % wages indexation
gammastar100=param(5);  % steady state growth rate of technology
gammamiu100=param(6);   % steady state growth rate of capital embodied technology
h=param(7);             % habit formation parameter
lambdapss=param(8);     % steady state of mark-up shock (pins down steady state of wages)
lambdaw=param(9);       % wage mark-up
Lss=param(10);          % steady state for leisure (the ss for leisure is pinned down by psiss. convenient to parameterize the model in terms of Lss)
pss100=param(11);       % steady state quarterly inflation (multiplied by 100)
Fbeta=param(12);        % weird parameter of SW that is a function of beta and the steady state quarterly real rate of interest (ensures that beta is less than 1)
gss=param(13);          % steady state of the public spending shock
niu=param(14);          % curvature of the utility function for leisure
xip=param(15);          % price stickiness
xiw=param(16);          % wage stickiness
chi=param(17);          % elasticity of the capital utilization cost function
S=param(18);            % investment adjustment cost
fp=param(19);           % reaction to annualized inflation
fy=param(20);           % reaction to annualized output growth 
rhoR=param(21);         % policy smoothing  
% AR coefficients of exogenous prorcesses 
rhoz=param(22); 
rhog=param(23); 
rhomiu=param(24); 
rholambdap=param(25); 
rhopsi=param(26); 
rhob=param(27); 
rhopitarg=param(28); 
rhomp=param(29); 
rhoupsilon=param(30); 
% Std of structural shocks 
sdstruct=param(31:40);
kappa       = -param(41);   % Scaling for spredd  
sdspreadme  =  param(42);   % Measurement error spread
% Std Policy Signals. Two vectors 
% stdsignals1 is based on the parameter estimates, expanded to cover
% NSIGNALS 
stdsignals=1e-7*ones(nsignals,1); 
stdsignals(1:4)=param(43:46); 

% stdsignals2 is for the counterfactuals
% This is added 2 SDX2 matrix below 
% Checked this is correct set-up for the filter in 
% zb_filterExpNomr.m 

% Alejandro Justiniano July 2011
if nargin < 2 || isfield(addin,'stdsignals')==0; 
    stdsignals2=1e-7*ones(nsignals,1); 
    disp('Setting STD all signals to 1e5'); 
    pause(0.5); 
else
    stdsignals2=addin.stdsignals;
    if length(stdsignals2)~=nsignals
        error('Wrong STD for signals')
    end
end 
% Parameters for Inflation Expectations
posCPI = param(47);
expInfLoading = param(48); 
end_nonp=48; 
% _________________________________________________________________________
% 3.2  Begin Parameters for Multiple Prices 
%__________________________________________________________________________
% 3.2.1  How to aggregate gdp deflator
flag_wadd1  =param(end_nonp+1); 
nserp       =param(end_nonp+2); 
% __________________________________________________________________
% 3.2.1 Declare positions for 
% Rows for 
% a) loadings 
% b) ar ID errors 
% c) constants 
% d) volatilities 
%
% *Note:* There are NSERP+1 loadings. First 2 are for the GDP Deflator.
% 1st of W(1)*Price consumption.2nd on W(2)*Prince Investment. 
% Remaining are for other prices on the Consumption Deflators 
row_load=(end_nonp+3):(end_nonp+nserp+3); 
row_ar  =row_load(end)*ones(1,nserp)+(1:nserp); 
row_std =row_ar(end)*ones(1,nserp)+(1:nserp);
row_c   =row_std(end)*ones(1,nserp)+(1:nserp); 

%__________________________________________________________________________
% 3.2.3 vector of factor loadings 
bet_gdpdef  =param(row_load(1:2)); 
bet_c       =param(row_load(3:end)); 

bet_CPI     =bet_c(posCPI-1); 
if abs(expInfLoading-999) < 0.01; 
    expInfLoading=bet_CPI;
end

%__________________________________________________________________________
% 3.2.4 vector of idiosyncratic errors 
rhoid=2*param(row_ar)-1; 
if any(rhoid < -0.995)~=0 || any(rhoid > 0.995)~= 0 
    disp('Explosive AR ID roots'); 
    GG=[];C=[];RR=[];eu=[0;0];SDX=[];ZZ=[];NY=[];NX=[];NETA=[];gev=[];ncof=[];
    addout=[]; 
    return
end

%__________________________________________________________________________
% 3.2.5 vector of constants 
cprices =param(row_c); 
if all( cprices ~= 0 ) 
    error('At least one of the constants must be set equal to 0')
end
constantCPI = cprices(posCPI);

%__________________________________________________________________________
% 3.2.6 vector of idiosyncratic errors
sdid             =param(row_std); 

%__________________________________________________________________________
% 3.2.7 constants for I deflator 
%% This is the end of the first parameter
cons_idef        =param(row_c(end)+1:ncof/2);
if length(cons_idef)> 1 
    error('Wrong Dimension of Investment Deflator') 
end 

%% paramS2 is the parameter for the second part of the sample
paramS2=param( (ncof/2)+1:end); 



%ncof=        62+4*(nserp-3);    % Number of coefficients 
%if length(param)~=ncof 
%    error('PARAM does not have NCOF entries') 
%end 

%==========================================================================
%% *PART 4. Compute Steady State (SS)* 
 
gammastar=gammastar100/100;
gammamiu=gammamiu100/100;
beta=100/(Fbeta+100);
rss=exp(gammastar)/beta-1;
rss100=rss*100;
gss=1/(1-gss);
Rkss=exp(gammastar+gammamiu)/beta-1+delta;
sss=1/(1+lambdapss);
wss=(sss*((1-alpha)^(1-alpha))/((alpha^(-alpha))*Rkss^alpha))^(1/(1-alpha));
kLss=(wss/Rkss)*alpha/(1-alpha);
FLss=(kLss^alpha-Rkss*kLss-wss);
yLss=kLss^alpha-FLss;

% Scale of the economy is fixed at 1, i.e. independent of LSS which now
% only enters observation equation for hours 
Lsscale=1;
kss=kLss*Lsscale;
iss=(1-(1-delta)*exp(-gammastar-gammamiu))*kss*exp(gammastar+gammamiu);
F=FLss*Lsscale;
yss=yLss*Lsscale;
css=yss/gss-iss;
    
%%==========================================================================

%==========================================================================
%% *PART 5. Declare System Matrices* 
% Model to be solved is written as 
% $$ A_{0}Y_{t+1} + A_{1}Y_{t}=A_{2}Y_{t-1} + PSI \eta_{t} $$  
A0  = zeros(NY,NY) ;
A1  = zeros(NY,NY) ;
A2  = zeros(NY,NY) ;
PSI = zeros(NY,NX) ;
%==========================================================================
expg=exp(gammastar);

%% *PART 6. Begin Filling Equation Rows* 

%% Eq 1, Output (y) from production function 
A1(y,y)=1;
A1(y,k)=-((yss+F)/yss)*alpha;
A1(y,L)=-((yss+F)/yss)*(1-alpha);

%% Eq 1*, Output (y*) in flex price & wage 
A1(ystar,ystar)=1;
A1(ystar,kstar)=-((yss+F)/yss)*alpha;
A1(ystar,Lstar)=-((yss+F)/yss)*(1-alpha);

%% Eq. 2, Effective Capital (K) 
A1(k,k)=1;
A1(k,u)=-1;
A1(k,z)=1;
A2(k,kbar)=1;
A1(k,miu)=alpha/(1-alpha)+1;

%% Eq 2*, Effective Capital (K*) in flex price & wage 
A1(kstar,kstar)=1;
A1(kstar,ustar)=-1;
A1(kstar,z)=1;
A2(kstar,kbarstar)=1;
A1(kstar,miu)=alpha/(1-alpha)+1;

%% Eq. 3 Hours (L), from cost minimization 
A1(L,Rk)=1;
A1(L,w)=-1;
A1(L,k)=1;
A1(L,L)=-1;

%% Eq. 3* Hours (L*) in flex price & wage 
A1(Lstar,Rkstar)=1;
A1(Lstar,wstar)=-1;
A1(Lstar,kstar)=1;
A1(Lstar,Lstar)=-1;

%% Eq. 4 Return to Capital (Rk) from utilization FOC
A1(Rk,Rk)=1;
A1(Rk,u)=-chi;

%% Eq. 4* Return to Capital (Rk*) in flex price & wage
A1(Rkstar,Rkstar)=1;
A1(Rkstar,ustar)=-chi;

%% Eq. 5, Real Wages (w) PC (with split of shocks)
slopewages=(1-beta*xiw)*(1-xiw)/((1+beta)*xiw*(1+niu*(1+1/lambdaw)));

dispaj('Slope in Wages:',slopewages); 
A1(w,w)=1;
A0(w,w)= -beta/(1+beta);
A0(w,p)= -beta/(1+beta);
A1(w,x)=slopewages;
A1(w,z)=(1+iotaw*beta-rhoz*beta)/(1+beta);
A1(w,miu) =(alpha/(1-alpha))*((1+iotaw*beta-rhomiu*beta)) *(1/(1+beta)); 

A1(w,p)=(1+iotaw*beta)/(1+beta);
A2(w,p)=iotaw/(1+beta);
A2(w,w)=1/(1+beta);

A2(w,z)=iotaw/(1+beta);
A2(w,miu)=(alpha*iotaw)/( (1+beta)*(1-alpha) ); 

A1(w,wmarkwn)= -1; 

%% Eq. 5*, Real Wages (w*) in flex price & wage 
% $$ x^{*}_{t} = 0 $$ i.e. wage gap is zero, so workers are always on their supply curve 
A1(wstar,xstar)=1;

%% Eq. 6, Wage Gap (x), (MRS-Real Wage)
A1(x,x)=1;
A1(x,w)=-1;
A1(x,lambda)=-1;
A1(x,b)=1/((1-rhob)*(expg-h*beta*rhob)*(expg-h)/(expg*h+expg^2+beta*h^2));
A1(x,L)=niu;
A1(x,psi)=1; 

%% Eq. 6*, Wage Gap (x*), (MRS-Real Wage) in flex price & wage 
% Equal zero from eq 5*
A1(xstar,xstar)=1;
A1(xstar,wstar)=-1;
A1(xstar,lambdastar)=-1;
A1(xstar,b)=1/((1-rhob)*(expg-h*beta*rhob)*(expg-h)/(expg*h+expg^2+beta*h^2));
A1(xstar,Lstar)=niu;
A1(xstar,psi)=1; 

%% Eq 7, Prices (p) PC 

dispaj('Slope in Prices:',((1-beta*xip)*(1-xip)/((1+iotap*beta)*xip)));
A1(p,p)=1;
A0(p,p)=-beta/(1+iotap*beta);
A2(p,p)=iotap/(1+iotap*beta);
A1(p,s)=-((1-beta*xip)*(1-xip)/((1+iotap*beta)*xip));
A1(p,lambdap)=-1;       

%% Eq 7*, Prices (p*) in flex price & wage
% $$ s^{*}_{t} = 0 $$ i.e. deviations of marginal cost from SS are zero, so
% firms are always on their supply curve
A1(Rstar,sstar)=1;

%% Eq 8, Marginal Cost (s)
A1(s,s)=1;
A1(s,Rk)=-alpha;
A1(s,w)=-(1-alpha);

%% Eq 8*, Marginal Cost (s*) in flex price & wage 
A1(sstar,sstar)=1;
A1(sstar,Rkstar)=-alpha;
A1(sstar,wstar)=-(1-alpha);

%% Eq 9, Lagrange Multiplier IBC (Lambda) 
A1(lambda,lambda)=1;
A1(lambda,R)=-1;
A0(lambda,lambda)=-1;
A0(lambda,p)=1;
A1(lambda,z)=rhoz;
A1(lambda,miu)=(rhomiu)*alpha/(1-alpha);

%% Eq 9*, Lagrange Multiplier (Lambda*) in flex price & wage
A1(lambdastar,lambdastar)=1;
A1(lambdastar,Rstar)=-1;
A0(lambdastar,lambdastar)=-1;
A1(lambdastar,z)=rhoz;
A1(lambdastar,miu)=(rhomiu)*alpha/(1-alpha);

%% Eq 10, Consumption (c)
A1(c,lambda)=(expg-h*beta)*(expg-h);
A1(c,b)=-(expg-h*beta*rhob)*(expg-h)/((1-rhob)*(expg-h*beta*rhob)*(expg-h)/(expg*h+expg^2+beta*h^2)); 
A1(c,z)=-(beta*h*expg*rhoz-h*expg);
A1(c,miu)=(-(beta*h*expg*rhomiu-h*expg))*alpha/(1-alpha);
A1(c,c)=(expg^2+beta*h^2);
A0(c,c)=-beta*h*expg;
A2(c,c)=expg*h;

%% Eq 10*, Consumption (c*) in flex price & wage
A1(cstar,lambdastar)=(expg-h*beta)*(expg-h);
A1(cstar,b)=-(expg-h*beta*rhob)*(expg-h)/((1-rhob)*(expg-h*beta*rhob)*(expg-h)/(expg*h+expg^2+beta*h^2));  
A1(cstar,z)=-(beta*h*expg*rhoz-h*expg);
A1(cstar,miu)=(-(beta*h*expg*rhomiu-h*expg))*alpha/(1-alpha);
A1(cstar,cstar)=(expg^2+beta*h^2);
A0(cstar,cstar)=-beta*h*expg;
A2(cstar,cstar)=expg*h;

%% Eq 11, Taylor Rule for Nominal Interest Rate (R)
A1(R,R)=1;
% 10.1 Systematic Component 
A2(R,R)=rhoR;
A1(R,dpma)  = -(1-rhoR)*fp; % Annualized Inflation 
A1(R,yobsan)= -(1-rhoR)*fy; % Annualized Observable Output 
% 10.2 Non-Systematic Component 
A1(R,mp)      = -1;         % Contemporaneous Innovation 
A1(R,sigR_0)  = -1;         % Signals 
A1(R,pitarg)  =(1-rhoR)*fp; % Inflation Drift

%% Eq. 12, Utilization (U) from market clearing condition 
A1(u,c)=css/yss;
A1(u,i)=iss/yss;
A1(u,y)=-1/gss;
A1(u,g)=1/gss;
A1(u,u)=kss*Rkss/yss;

%% Eq. 12*, Utilization (U) in flex price and wage 
A1(ustar,cstar)=css/yss;
A1(ustar,istar)=iss/yss;
A1(ustar,ystar)=-1/gss;
A1(ustar,g)=1/gss;
A1(ustar,ustar)=kss*Rkss/yss;

%% Eq. 13, Shadow Value of Capital (phi)
A1(phi,phi)=1;
A0(phi,phi)=-beta*exp(-gammastar-gammamiu)*(1-delta);
A1(phi,z)=rhoz;
A1(phi,miu)=(rhomiu)*(alpha/(1-alpha)+1);
A0(phi,lambda)=-(1-beta*exp(-gammastar-gammamiu)*(1-delta));
A0(phi,Rk)=-(1-beta*exp(-gammastar-gammamiu)*(1-delta));

%% Eq. 13*, Shadow Value of Capital (phi*) in flex price and wage
A1(phistar,phistar)=1;
A0(phistar,phistar)=-beta*exp(-gammastar-gammamiu)*(1-delta);
A1(phistar,z)=rhoz;
A1(phistar,miu)=(rhomiu)*(alpha/(1-alpha)+1);
A0(phistar,lambdastar)=-(1-beta*exp(-gammastar-gammamiu)*(1-delta));
A0(phistar,Rkstar)=-(1-beta*exp(-gammastar-gammamiu)*(1-delta));

%% Eq. 14, Investment (i, in consumption units) 
expgmiu=exp(gammastar+gammamiu);
A1(i,lambda)=1/(S*expgmiu^2);
A1(i,phi)=-1/(S*expgmiu^2);
A1(i,upsilon)=-1/(S*expgmiu^2);
A1(i,miu)=((1-beta*rhomiu))*(alpha/(1-alpha)+1); %no normalization
A1(i,i)=(1+beta);
A1(i,z)=(1-beta*rhoz);
A0(i,i)=-beta;
A2(i,i)=1;

%% Eq. 14*, Investment (i*, in consumption units) in flex price and wage 
A1(istar,lambdastar)=1/(S*expgmiu^2);
A1(istar,phistar)=-1/(S*expgmiu^2);
A1(istar,upsilon)=-1/(S*expgmiu^2);
A1(istar,miu)=((1-beta*rhomiu))*(alpha/(1-alpha)+1); %no normalization
A1(istar,istar)=(1+beta);
A1(istar,z)=(1-beta*rhoz);
A0(istar,istar)=-beta;
A2(istar,istar)=1;

%% Eq 15, Capital accumulation (kbar)
% Note: Recall, effective capital above is KBAR*UTILIZATION 
A1(kbar,kbar)=1;
A1(kbar,i)=-(1-(1-delta)*exp(-gammastar-gammamiu));
A1(kbar,upsilon)=-(1-(1-delta)*exp(-gammastar-gammamiu));
A2(kbar,kbar)=(1-delta)*exp(-gammastar-gammamiu);
A1(kbar,z)=(1-delta)*exp(-gammastar-gammamiu);
A1(kbar,miu)=((1-delta)*exp(-gammastar-gammamiu))*(alpha/(1-alpha)+1);

%% Eq 15*, Capital accumulation (kbar*) in flex price and wage 
A1(kbarstar,kbarstar)=1;
A1(kbarstar,istar)=-(1-(1-delta)*exp(-gammastar-gammamiu));
A1(kbarstar,upsilon)=-(1-(1-delta)*exp(-gammastar-gammamiu));
A2(kbarstar,kbarstar)=(1-delta)*exp(-gammastar-gammamiu);
A1(kbarstar,z)=(1-delta)*exp(-gammastar-gammamiu);
A1(kbarstar,miu)=((1-delta)*exp(-gammastar-gammamiu))*(alpha/(1-alpha)+1);

%% Eq 16, MP Exogenous Policy Innovation (unexpected) 
A1(mp,mp)=1; A2(mp,mp)=rhomp; PSI(mp,Rs)=1; 

%% Eq 17, Exogenous Growth Rate in Neutral Technoloy 
A1(z,z)=1; A2(z,z)=rhoz; PSI(z,zs)=1;

%% Eq 18, Exogenous Process in G Spending 
A1(g,g)=1; A2(g,g)=rhog; PSI(g,gs)=1;

%% Eq 19, Exogenous Process in Investment Specific Tech or Relative Price C
% to I 
A1(miu,miu)=1; A2(miu,miu)=rhomiu; PSI(miu,mius)=1;

%% Eq 20, Exogenous Process in Price Markup 
A1(lambdap,lambdap)=1; A2(lambdap,lambdap)=rholambdap;PSI(lambdap,lambdaps)=1;

%% Eq 21, Exogenous Process in Labor Disutility (AR1) 
A1(psi,psi)=1; A2(psi,psi)=rhopsi; PSI(psi,psis)=1;

%% Eq 22, Exogenous Process in Intertemporal Preference
A1(b,b)=1; A2(b,b)=rhob; PSI(b,bs)=1;

%% Eq 23, Exogenous Process in MEI 
A1(upsilon,upsilon)=1; A2(upsilon,upsilon)=rhoupsilon; PSI(upsilon,upsilons)=1;

%% Eq 24, Exogenous Process in Inflation Drift 
A1(pitarg,pitarg)=1; A2(pitarg,pitarg)=rhopitarg; PSI(pitarg,shpitarg)=1;

%% Eq 25, Exogenous Process in Wage Markup (white noise) 
A1(wmarkwn,wmarkwn)=1; PSI(wmarkwn,wmarkwnsh)=1; 

%% Eq. 26, Observable GDP Nominal (add P back)
A1(ynomobs,ynomobs) =  1;
A1(ynomobs,gdp)     = -1; 
A1(ynomobs,dagtrend)= -1; 
% Add inflation, since nominal 
A1(ynomobs,p)       = -1; 

A2(ynomobs,gdp)     = -1; 

%% Eq. 27, Observable Consumption Nominal (add P back)
A1(cnomobs,cnomobs) =1;
A1(cnomobs,c)       = -1; 
A1(cnomobs,dagtrend)= -1; 
% Add inflation, since nominal 
A1(cnomobs,p)= -1; 

A2(cnomobs,c)= -1; 

%% Eq. 28 Observable Investment Nominal (add P back)
% In consumption units
A1(inomobs,inomobs)   =1;
A1(inomobs,i)         = -1; 
A1(inomobs,dagtrend)  = -1; 
% Add inflation, since nominal 
A1(inomobs,p)= -1; 

A2(inomobs,i)= -1; 

%% Eq. 29 Observable Hourly Wages  Nominal (add P back)
A1(wnomobs,wnomobs) =1;
A1(wnomobs,w)       = -1; 
A1(wnomobs,dagtrend)= -1; 
% Add inflation, since nominal 
A1(wnomobs,p)= -1; 

A2(wnomobs,w)= -1; 

%% Eq. 30 Observable Real GDP Growth (in C Units)
% Corrected observation equation for utilization 
% $$ yrobs_{t} = gdp_{t}-gdp_{t-1} + agtrend_{t} $$ 
% where AGTREND is the aggregate trend 
A1(yrobs,yrobs)=1; 
A1(yrobs,gdp)= -1;
A1(yrobs,dagtrend)= -1; 
A2(yrobs,gdp)= -1; 

%% Eq. 31 GDP (detrended), Real Detrended Output corrected for Utilization 
A1(gdp,gdp)= -1; 
A1(gdp,y)=1; 
A1(gdp,u)= -kss*Rkss/yss; 

%% Eq. 31* GDP* in flex price and wage
A1(gdpstar,gdpstar)= -1; 
A1(gdpstar,ystar)=1; 
A1(gdpstar,ustar)= -kss*Rkss/yss; 

%% Eq. 32 Growth Rate in Investment Deflator 
% $$ \Delta P^{I}_{t} = \Delta P^{C}_{t} - \mu_{t} $$
A1(dpinv,dpinv)= -1; 
A1(dpinv,p)    =  1; 
A1(dpinv,miu)  = -1; 

%% Eq. 33 GDP Deflator 
% $$ P^{GDP}_{t} = \beta_{1} P^{C}_{t} + \beta_{2} P^{I}_{t} $$ 
% 
% Recall there are 2 ways to define the weights 
%
% *First:* $$ {\frac{C}{C+I} , \frac{I}{C+I} } $$ 
% 
% *Second:* $$ {\frac{C}{Y}  , \frac{I}{Y}   } $$
%
% Note that in second case weights will not add up to one, which is not
% desirable, but opens up the door for adding a 3 deflator, i.e. price of G
% and NX. 
bet_gdpdef=bet_gdpdef(:); 
weights=zeros(1,2); 
if flag_wadd1==1
    weights(1)=css/(css+iss);
    weights(2)=iss/(css+iss);
else
    weights(1)=css/yss;
    weights(2)=iss/yss;
end
A1(gdpdef,gdpdef) = -1; 
A1(gdpdef,p)      = weights(1)*bet_gdpdef(1); 
A1(gdpdef,dpinv)  = weights(2)*bet_gdpdef(2); 


%% Eq. 34, Nominal C deflated by Model-based Pc deflator 
A1(crobs,crobs)  =-1; 
A1(crobs,cnomobs)=1; 
A1(crobs,p      )= -1; 
% This expression is easier than 
% A1(crobs,crobs)=1; 
% A1(crobs,c)= -1;
% A1(crobs,z)= -1; 
% A1(crobs,miu)= -alpha/(1-alpha); 
% A2(crobs,c)= -1; 

%% Eq. 35, Nominal I deflated by Model-based Pc deflator 
A1(irobs,irobs)  =-1; 
A1(irobs,inomobs)=1; 
A1(irobs,p      )= -1; 

%% Eq. 36, Nominal I deflated by Model-based Pc deflator 
A1(wrobs,wrobs)  =-1; 
A1(wrobs,wnomobs)=1; 
A1(wrobs,p      )= -1; 

%% Eq. 37, Output gap $$ GDP_{t} - GDP^{*}_{t} $$ 
A1(ygap,ygap)     =-1; 
A1(ygap,gdp)      =1; 
A1(ygap,gdpstar)  = -1; 
 
%% Eq 38 EXREAL1: Expected Real Rate next period $$ E[ R_{t} - \pi_{t+1} | t ] $$
A1(exreal1,exreal1) = -1; 
A1(exreal1,R)       =1; 
A0(exreal1,p)       = -1; 

%% Eq. 39 Lag 1 Observable Real GDP Growth (in C Units)
A1(yrobs_1,yrobs_1)=1; A2(yrobs_1,yrobs)=1;

%% Eq. 40 Lag 2 Observable Real GDP Growth (in C Units)
A1(yrobs_2,yrobs_2)=1; A2(yrobs_2,yrobs_1)=1;

%% Eq. 41 Annualized Observable Real GDP Growth (in C Units)
% $$ \Delta y^{A}_{t}= \sum_{j=0}^{3}yrobs_{t-j} $$
A1(yobsan,yobsan)= -1;

A1(yobsan,yrobs)  =1; 
A1(yobsan,yrobs_1)=1; 
A1(yobsan,yrobs_2)=1;
A2(yobsan,yrobs_2)= -1;

%% Eq. 42 Lag 1 Inflation (C)
A1(p_1,p_1)=1; A2(p_1,p)=1;
%% Eq. 43 Lag 2 Inflation (C)
A1(p_2,p_2)=1; A2(p_2,p_1)=1;
%% Eq. 44 Annualized Inflation (C) 
% $$ \Delta P^{A,C}_{t}= \frac{ \sum_{j=0}^{3} \Delta P^{C}_{t-j} }{4} $$
A1(dpma,dpma)= -1;
A1(dpma,p)   = 1/4;
A1(dpma,p_1) = 1/4;
A1(dpma,p_2) = 1/4;
A2(dpma,p_2) = -1/4;
%% Eq. 45, Growth Rate in Aggregate Trend 
% $$ \Gamma_{t} = Z_{t} + \frac{\alpha}{1-\alpha} \mu_{t} $$ 
A1(dagtrend,dagtrend)= -1; 
A1(dagtrend,z)=1; 
A1(dagtrend,miu)=alpha/(1-alpha); 

%% Eq. 46-58 Expectations of Future Interest Rates 
% $$ E[ R_{t+j} | t ] j=1,..,4 $$ In that order 
A1(R_lead1,R_lead1)=1; A0(R_lead1,R)= -1; 
A1(R_lead2,R_lead2)=1; A0(R_lead2,R_lead1)= -1;
A1(R_lead3,R_lead3)=1; A0(R_lead3,R_lead2)= -1;
A1(R_lead4,R_lead4)=1; A0(R_lead4,R_lead3)= -1;
A1(R_lead5,R_lead5)=1; A0(R_lead5,R_lead4)= -1;
A1(R_lead6,R_lead6)=1; A0(R_lead6,R_lead5)= -1;
A1(R_lead7,R_lead7)=1; A0(R_lead7,R_lead6)= -1;
A1(R_lead8,R_lead8)=1; A0(R_lead8,R_lead7)= -1;
A1(R_lead9,R_lead9)=1; A0(R_lead9,R_lead8)= -1;
A1(R_lead10,R_lead10)=1; A0(R_lead10,R_lead9)= -1;
A1(R_lead11,R_lead11)=1; A0(R_lead11,R_lead10)= -1;
A1(R_lead12,R_lead12)=1; A0(R_lead12,R_lead11)= -1;


%% Eq. 45-49 Monetary Policy Signals

%% Eq. $$ \zeta^{12}_{t}       = \epsilon^{12}_{t} $$ 
A1(sigR_12,sigR_12)=1;PSI(sigR_4,Rs_lag12)=1; 

%% Eq. $$ \zeta^{11}_{t}=\zeta^{12}_{t-1} + \epsilon^{11}_{t} $$ 
A1(sigR_11,sigR_11)=1;A2(sigR_11,sigR_12)=1;PSI(sigR_11,Rs_lag11)=1; 

%% Eq. $$ \zeta^{10}_{t}=\zeta^{11}_{t-1} + \epsilon^{10}_{t} $$ 
A1(sigR_10,sigR_10)=1;A2(sigR_10,sigR_11)=1;PSI(sigR_10,Rs_lag10)=1; 

%% Eq. $$ \zeta^{9}_{t}=\zeta^{10}_{t-1} + \epsilon^{9}_{t} $$ 
A1(sigR_9,sigR_9)=1;A2(sigR_9,sigR_10)=1;PSI(sigR_9,Rs_lag9)=1; 

%% Eq. $$ \zeta^{8}_{t}=\zeta^{9}_{t-1} + \epsilon^{8}_{t} $$ 
A1(sigR_8,sigR_8)=1;A2(sigR_8,sigR_9)=1;PSI(sigR_8,Rs_lag8)=1; 

%% Eq. $$ \zeta^{7}_{t}=\zeta^{8}_{t-1} + \epsilon^{7}_{t} $$ 
A1(sigR_7,sigR_7)=1;A2(sigR_7,sigR_8)=1;PSI(sigR_7,Rs_lag7)=1; 

%% Eq. $$ \zeta^{6}_{t}=\zeta^{7}_{t-1} + \epsilon^{6}_{t} $$ 
A1(sigR_6,sigR_6)=1;A2(sigR_6,sigR_7)=1;PSI(sigR_6,Rs_lag6)=1; 

%% Eq. $$ \zeta^{5}_{t}=\zeta^{6}_{t-1} + \epsilon^{5}_{t} $$ 
A1(sigR_5,sigR_5)=1;A2(sigR_5,sigR_6)=1;PSI(sigR_5,Rs_lag5)=1; 

%% Eq. $$ \zeta^{4}_{t}=\zeta^{5}_{t-1} + \epsilon^{4}_{t} $$ 
A1(sigR_4,sigR_4)=1;A2(sigR_4,sigR_5)=1;PSI(sigR_4,Rs_lag4)=1; 

%% Eq. $$ \zeta^{3}_{t} = \zeta^{4}_{t-1} + \epsilon^{3}_{t} $$ 
A1(sigR_3,sigR_3)=1; A2(sigR_3,sigR_4)=1; PSI(sigR_3,Rs_lag3)=1; 

%% Eq. $$ \zeta^{2}_{t} = \zeta^{3}_{t-1} + \epsilon^{2}_{t} $$ 
A1(sigR_2,sigR_2)=1; A2(sigR_2,sigR_3)=1; PSI(sigR_2,Rs_lag2)=1; 

%% Eq. $$ \zeta^{1}_{t} = \zeta^{2}_{t-1} + \epsilon^{1}_{t} $$ 
A1(sigR_1,sigR_1)=1; A2(sigR_1,sigR_2)=1; PSI(sigR_1,Rs_lag1)=1; 

%% Eq. $$ \zeta^{0}_{t} = \zeta^{1}_{t-1} $$, which implies that the history of
% policy signals regarding the current quarter is given by 
%
% $$ \zeta^{0}_{t} = \sum_{j=1}^{12} \epsilon^{j}_{t-j}  $$
%
% *Note:* The contemporaneous policy innovation appears 
% directly in the definition of the Taylor rule, in case it is AR(1). An
% alternative would be to make the current composite signal +
% contemporaneous innovation persistent. 
A1(sigR_0,sigR_0)=1;A2(sigR_0,sigR_1)=1;  

%% Inflation Expectations Equation
pe_vec = pe_vec(:);
A1(pe_vec,pe_vec)=eye(length(pe_vec)); 
A0(pe_vec,[p;pe_vec(1:end-1)])= -eye(length(pe_vec)); 

A1(expectedPi40,expectedPi40)= -1; 
A1(expectedPi40,pe_vec      )= (expInfLoading/length(pe_vec))*ones(1,length(pe_vec));  

%% Sum of Real Rates Equation
A1(realRate12q,realRate12q)= -1; 
A1(realRate12q,R          )=(1/12)*1; 
A1(realRate12q,Rlead_vec(1:11) )=(1/12)*ones(11,1); 
A1(realRate12q,pe_vec(1:12) )= (1/12)*-ones(12,1);

%% *PART 7: Solve Model using AMASOLVE*
[GG,~,impact,~,~,~,gev,eu]=amasolve(A0,A1,A2,zeros(NY,1),PSI) ;
% EU =[1;1] if unique and stable RE solution 
if ~isequal(eu,[1;1]); 
    ZZ=[]; C=[]; SDX=[];
    return 
end 

%% *PART 8: Addition of Measurement Equations for Prices*

%% 8.1 Extend model solution by 2*NSERP rows 
% Number of rows to add to the solution 
nadd=nserp*2; 
% Part 1: Position of prices 
posp =(NY+1):(NY+nserp); 
% Part 2: Position of ID errors 
posid=(posp(end)+1):(NY+2*nserp); 
GG=blkdiag(GG,zeros(nadd)); 
% Matrix of impact coefficients 
RR=zeros(NY+nadd,NX+nserp); 
RR(1:NY,1:NX)=impact; 
clear impact; 
%% 8.2 Idiosyncratic errors 
% Which evolve as $$ \nu^{i}_{t} = \rho  \nu^{i}_{t-1} + \varsigma^{i}_{t} $$
% as independent AR(1) processes 
GG(posid,posid   )=diag(rhoid); 
RR(posid,NX+1:end)=eye(nserp); 

%% 8.3 General Structure of Observable Price Equation 
% 
% For price $$ P^{i,C}_{t} = \beta^{i} P^{C}_{t} + \nu^{i}_{t} $$ written in terms
% of lagged variables as 
% $$ P^{i,C}_{t} = \beta^{i} G_{P^{C}} S_{t-1} + \beta^{i} R_{P^{C}}
% \eta_{t} + \rho \nu^{i}_{t-1} + \varsigma^{i}_{t} $$ ,
% where $$ \beta^{i} $$ is the series specific factor loading, and, $$ G_{P^{C}}  R_{P^{C}} $$ correspond to the rows of the solution
% matrix for the model consumption price 

%% 8.3.a First row is GDP 
GG(posp(1),1:NY)=GG(gdpdef,1:NY);
RR(posp(1),1:NX)=RR(gdpdef,1:NX); 

%% 8.3.b Remaining rows are multiple indicators for consumption prices 
% Recall Structural part of observation equation $$ \beta^{i} G_{P^{C}} S_{t-1} +
% \beta^{i} R_{P^{C}} \eta_{t}  $$ 
GG(posp(2:end),:)       =repmat( GG(p,:)    ,[nserp-1 1] ); 
RR(posp(2:end),1:NX)    =repmat( RR(p,1:NX) ,[nserp-1 1] );
% Multiply by Factor Loadings  
bet_c=bet_c(:); 
GG(posp(2:end),1:NY)=(repmat(bet_c,[1 NY])).*GG(posp(2:end),1:NY); 
RR(posp(2:end),1:NX)=(repmat(bet_c,[1 NX])).*RR(posp(2:end),1:NX);
%% 8.3.c Add idiosyncratic errors 
GG(posp,posid)   =diag(rhoid); 
RR(posp,NX+1:end)=eye(nserp); 

% Update model dimensions after adding prices 
%impact=RR; 
NY=size(GG,1); 
NX=size(RR,2); 

%% *Part 9: Addition of Measurement Equations for Spread*
% 
%% 9.1 Extend model solution by 1 row 
GG    =blkdiag(GG,zeros(2));
RR   =[RR zeros(NY,1);zeros(2,NX+1)];

%% 9.2 Fill in last row for the spread 
% 
% $$ sp_{t} = \kappa \upsilon_{t} + \varsigma^{sp}_{t} $$ 
posspread=NY+1;
GG(posspread,:)   =kappa*GG(upsilon,:); 
RR(posspread,:)   =kappa*RR(upsilon,:); 
RR(posspread,NX+1)=1; 

%% 9.3 Fill in spread measurement error
% 
RR(posspread+1,NX+1)=1; 

% Update model dimensions after adding spread
NY=NY+2;
NX=NX+1;


%% *Part 10: Declare Real Observable {Y,C,I,W,SumR}, Accounting for Measurement Errors* 
nadd2=4; 
GG       =blkdiag(GG,zeros(nadd2)); 
RR       =[RR;zeros(nadd2,size(RR,2))]; 
posGDPreal  =NY+1; 
posCreal    =NY+2; 
posIreal    =NY+3; 
posWreal    =NY+4; 
%posSumR    =NY+5;

%% Eq. Real GDP: Nominal Deflated by GDP deflator 
GG(posGDPreal,:)     =GG(ynomobs,: )-GG( posp(1),:); 
RR(posGDPreal,:)     =RR(ynomobs,: )-RR( posp(1),:); 

%% Eq. Real Consumption: Nominal Deflated by C deflator 
GG(posCreal,:) =GG(cnomobs,: )-GG( posp(3),:); 
RR(posCreal,:) =RR(cnomobs,: )-RR( posp(3),:); 

%% Eq. Real Investment: Nominal Deflated by I deflator 
GG(posIreal,:) =GG(inomobs,: )-GG(dpinv,:); 
RR(posIreal,:) =RR(inomobs,:) -RR(dpinv,:); 

%% Eq. Real Wages: Nominal Deflated by C deflator 
GG(posWreal,:) =GG(wnomobs,: )-GG(posp(3),:); 
RR(posWreal,:) =RR(wnomobs,: )-RR(posp(3),:); 

%% Sum of Expected Real Rates over 3 Year Horizon
% GG(posSumR,:) = GG(R_lead1,:)+GG(R_lead2,:)+GG(R_lead3,:)+GG(R_lead4,:)+GG(R_lead5,:)+GG(R_lead6,:)+...
%                 GG(R_lead7,:)+GG(R_lead8,:)+GG(R_lead9,:)+GG(R_lead10,:)+GG(R_lead11,:)+GG(R_lead12,:)-...
%                 GG(pe_vec(2),:)-GG(pe_vec(3),:)-GG(pe_vec(4),:)-GG(pe_vec(5),:)-GG(pe_vec(6),:)-GG(pe_vec(7),:)-GG(pe_vec(8),:)-...
%                 GG(pe_vec(9),:)-GG(pe_vec(10),:)-GG(pe_vec(11),:)-GG(pe_vec(12),:)-GG(pe_vec(13),:);  
NY=size(GG,1); 

%% *PART 11: Declare Cholesky of Variance Covariance Matrix* 
% Recall $$ V = SDX'SDX $$ and the order of the shocks is 
% 1. Structural 2. Policy Signals 3. Idiosyncratic Prices 4. Measurement
% Error Spread 
SDX=diag([sdstruct(:);stdsignals;sdid(:);sdspreadme]);
if size(SDX,1)~=NX;error('Dimensions of errors and SDX do not match');end

%% *PART 12: Declare ZZ, Matrix of Pointers to Observables (Estimation 1)*
% Ordering 
% 1) Nominal GDP Growth. 
% 2) Nominal C   Growth.
% 3) Nominal I   Growth.
% 4) Hours. 
% 5) Nominal W   Growth. 
% 6) Feds Fund Rate. 
% 7) Investment Deflator. 
% 8) GDP Deflator. 
% 9) through NSERP-1) Indicators of C Prices. 
% 10) Spread
% Last is 10 yr Avg. Expected Inflation
dimZ=7+nserp+1+1; 
ZZ=zeros(dimZ,NY); 
ZZ(1,ynomobs)=1; 
ZZ(2,cnomobs)=1; 
ZZ(3,inomobs)=1; 
ZZ(4,L)=1; 
ZZ(5,wnomobs)=1; 
ZZ(6,R)=1;
ZZ(7,dpinv)=1;
ZZ(8:dimZ-2,posp)=eye(nserp);
ZZ(dimZ-1,posspread)=1; 
ZZ(dimZ  ,expectedPi40)=1;  

%% *PART 13: Declare C, Matrix of Constants*
C=zeros(NY,1); 
C(ynomobs)=gammastar100+pss100;
C(cnomobs)=gammastar100+pss100;
C(inomobs)=gammastar100+pss100;
C(L)=Lss;
C(wnomobs)=gammastar100+pss100;
C(R)    =pss100+rss100;
C(dpinv)=(pss100-gammamiu100)+cons_idef; 
C(posp)=pss100*ones(nserp,1)+cprices(:); 
C(expectedPi40)=pss100+constantCPI;
%ncof=ncof+1; 

%% *PART 14: Declare Additional Z2 Matrix (ADDOUT.Z2)*  
% Includes automatic rows for signals turned on 
PosActiveSig=find( stdsignals2 > 0.02 ); 
nActiveSig  =length( PosActiveSig );
dispaj('Number of Automatic signals',nActiveSig); 
if nActiveSig > 1
    addout.Z2=zeros(nActiveSig,NY);
    addout.Z2(1:nActiveSig, Rlead_vec(1:nActiveSig) )=diag( ones(nActiveSig,1) );
else
    addout.Z2=[];
end
%% *PART 15: Declare Additional addout.C2 Matrix (ADDOUT.addout.C2)*  
% Used to automatically fill in constants when forecasting 
addout.C2=C; 
addout.C2(p)=pss100; 
addout.C2(pitarg)=pss100; 
addout.C2(yrobs) =gammastar100; 
addout.C2(gdpdef)=pss100; 
addout.C2(crobs) =gammastar100; 
addout.C2(irobs) =gammastar100; 
addout.C2(wrobs) =gammastar100; 
addout.C2(exreal1)=rss100; 
addout.C2(yrobs_1) =gammastar100; 
addout.C2(yrobs_2) =gammastar100; 
addout.C2(yobsan) =gammastar100;
addout.C2([p_1;p_2;dpma])=pss100; 
addout.C2(dagtrend)=gammastar100; 
addout.C2(R_lead1:R_lead12)=pss100+rss100; 
addout.C2(posGDPreal)=gammastar100+pss100-C(posp(1)); 
addout.C2(posCreal  )=gammastar100+pss100-C(posp(2)); 
addout.C2(posIreal  )=gammastar100+C(dpinv); 
addout.C2(posWreal  )=gammastar100+pss100-C(posp(2)); 

addout.SDX2=diag([sdstruct(:);stdsignals2;sdid(:);sdspreadme]);

%% *PART 16: Declare Names of States and Shocks*

%% 16.1 STANAMES: Cell with Names of ALL states 

% Last Modified 
% AJ July 26 2011 
stanames1={'yhat','khat','hours','Return K','what','x','dp','Marginal Cost','lambdahat','chat',...
    'nomr','Utilization ','phihat','ihat','kbar','mp','z','g','miu','lambdap',...
    'psi','b','upsilon ','pitarg','wmarkwn','GDPm nominal','Cm nominal','Im nominal','Wm nominal','GDPm real',...
    'GDPm level','I deflator','GDPm deflator','Cm real','Im real','Wm real','Output gap','Real rate','GDPm lag1','GDPm lag2',...
    'GDP annualized','Inflation lag1','Inflation lag2','Inflation annualized','Aggregate trend','ExFFR1','ExFFR2','ExFFR3','ExFFR4',...
    'ExFFR5','ExFFR6','ExFFR7','ExFFR8','ExFFR9','ExFFR10','ExFFR11','ExFFR12','Signal0','Signal1','Signal2',...
    'Signal3','Signal4','Signal5','Signal6','Signal7','Signal8','Signal9','Signal10','Signal11','Signal12',...
    'ExpInf1','ExpInf2','ExpInf3','ExpInf4','ExpInf5','ExpInf6','ExpInf7','ExpInf8','ExpInf9','ExpInf10',...
    'ExpInf11','ExpInf12','ExpInf13','ExpInf14','ExpInf15','ExpInf16','ExpInf17','ExpInf18','ExpInf19','ExpInf20',...
    'ExpInf21','ExpInf22','ExpInf23','ExpInf24','ExpInf25','ExpInf26','ExpInf27','ExpInf28','ExpInf29','ExpInf30',...
    'ExpInf31','ExpInf32','ExpInf33','ExpInf34','ExpInf35','ExpInf36','ExpInf37','ExpInf38','ExpInf39','ExpInf40','AvgExpInf10y','RealRate12q',...
    'ystar','kstar','hstar','rkstar','wstar','xstar','rstar','mgcstar','lambdastar','cstar',...
    'ustar','phistar','istar','kbarstar','GDP star','GDP deflator'}; 
if nserp==3 
    stanames2={'C deflator','PCE Core','GDP def Id','C def Id','PCE Core Id'}; 
elseif nserp==4 
    stanames2={'C deflator','PCE Core','CPI Core','GDP def Id','C def Id','PCE Core Id','CPI Core Id'}; 
else 
    disp('Nprices > 4 please add names to STANAMES') 
end 
stanames3={'Spread','Spread Me','RealY','RealC','RealI','RealW'}; 
addout.stanames=[stanames1(:);stanames2(:);stanames3(:)]; 

%% 16.2 STSHOCKS: Cell with Names of Shocks in STANAMES 
% *Note:* Technically not true for the policy signals, for which looking
% only at the cummulative signals. Hence isolate policy innovations rather
% than policy states.
addout.stshocks={'mp','z','g','miu','lambdap','psi','b','upsilon ','pitarg','wmarkwn',...
    'Signal1','Signal2','Signal3','Signal4','Signal5','Signal6','Signal7','Signal8','Signal9','Signal10','Signal11','Signal12',...
    'GDP def Id','C def Id','PCE Core Id'};
if nserp==4
    addout.stshocks=[addout.stshocks,'CPI Core Id'];
end
addout.stshocks=[addout.stshocks,'Spread Me'];
addout.stshocks=addout.stshocks(:);

%% 16.3 SHONAMES: Cell with Names of Shock Innovations 
addout.shonames={'MP','Neutral','G+NX','ISTS','Price Markup','Labor Disutility',...
    'Discount','MEI','Inflation Drift','Wage Markup',...
    'MP Signal 1','MP Signal 2','MP Signal 3','MP Signal 4','MP Signal 5','MP Signal 6','MP Signal 7',...
    'MP Signal 8','MP Signal 9','MP Signal 10','MP Signal 11','MP Signal 12'}; 
addout.shonames=[addout.shonames(:);addout.stshocks(23:end)]; 

%% 16.4 Quick check dimensions 
if length(addout.stanames)~=size(GG,1)
    error('Mistmatch STANAMES and NY');
end
if length(addout.stshocks)~=size(SDX,1)
    error('Mistmatch STSHOCKS and NX');
end
if length(addout.shonames)~=size(SDX,1)
    error('Mistmatch SHONAMES and NX');
end
